#ifndef COMB_CARTESIAN_PRODUCT_HPP__
#define COMB_CARTESIAN_PRODUCT_HPP__

#include <vector>
#include <cmath>
#include <iostream>

/** Generate a vector of all ordered combinations (with replacement) of the members of elements, cnt times.
 *
 * called with ( {1,2} 2 ) this will generate the list {1,1},{1,2},{2,1},{2,2}
 *
 * \param elements An array of elements to be used.
 * \param cnt The length of the sequences to be generated.
 */
template<typename T>
std::vector<std::vector<T> > cartesian_product( std::vector<T> &elements, unsigned int cnt)
{
    if(cnt == 0) 
        error("cnt must be > 0 in cartesian_product");

    unsigned int i, j, N, total;
    
    N = elements.size();
    
    // calculate the number of elements in the cartesian product.
    total = (unsigned int)pow((double)N, cnt);

    std::vector<std::vector<T> > cprod;
    cprod.resize(total);

    // fill the vector one column at a time.  The frequency at which the
    // 'odometer' moves depends on the position of the colunm we are processing.
    unsigned int speed = total;
    unsigned int len = elements.size();
    // for each column of the output
    for(i = 0; i < cnt; i++)
    {
        speed /= N;
        // push one element onto the jth block. 
        for(j = 0; j < total ; j++)
            cprod[j].push_back(elements[j / speed % len]);
    }
    return cprod;
}

/** Calculate the cartesian product of several lists of items. 
 *
 * with {{1,2}, {3,4}} this will produce {1,3},{1,4},{2,3},{2,4}
 *
 * \param elements List of lists of elements to be expanded.
 */
template<typename T>
std::vector<std::vector<T> > cartesian_product( std::vector<std::vector<T> > &elements)
{
    unsigned int i,j,total;
    
    // calculate the number of elements in the cartesian product.
    total = 1;

    for(i = 0; i < elements.size(); i++)
        total *= elements[i].size();
    
    std::vector<std::vector<T> > cprod;
    cprod.resize(total);
    
    int speed = total;

    for(i = 0; i < elements.size(); i++)
    {
        int len = elements[i].size();
        speed /= len;

        for(j = 0; j < total ; j++)
            cprod[j].push_back(elements[i][j / speed % len]);
    }
    return cprod;
}

#endif
